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Abstract 

Si 

CN | A Cauchy-characteristic initial value problem for the Einstein-Klein- 

\ Gordon system with spherical symmetry is presented. Initial data are spec- 

ified on the union of a space-like and null hypersurface. The development 
of the data is obtained with the combination of a constrained Cauchy evo- 
lution in the interior domain and a characteristic evolution in the exterior, 
asymptotically flat region. The matching interface between the space-like 
\Q \ and characteristic foliations is constructed by imposing continuity conditions 

on metric, extrinsic curvature and scalar field variables, ensuring smoothness 
across the matching surface. The accuracy of the method is established for 
i ■ all ranges of M/R, most notably, with a detailed comparison of invariant ob- 

(3JT), servables against reference solutions obtained with a calibrated, global, null 

j> | algorithm. 

•i-H . 

X: 

o3 ■ I. INTRODUCTION 

The correct physical formulation of any asymptotically flat, radiative Cauchy problem 
requires boundary conditions at spatial infinity. These conditions ensure not only that total 
energy and the energy loss by radiation are both finite, but they are also responsible for 
the proper 1/r asymptotic falloff of the radiation fields. However, when treating radiative 
systems computationally, an outer boundary must be established artificially at some large 
but finite distance in the wave zone, i.e. many wavelengths from the source. Imposing an 
accurate radiation boundary condition at a finite distance is a difficult task even in the 
case of simple radiative systems evolving on a fixed geometric background. The problem is 
exacerbated when dealing with the Einstein equations. 

In recent years, the characteristic initial value problem (cIVP) formulation of the Ein- 
stein equations in non-spherically symmetric configurations [|18| , |15| has been explored, pro- 
viding possible alternatives to the practical and theoretical problems introduced by the outer 
boundary conditions of the Cauchy initial value problem (CIVP). Based on the concept of 
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combined Cauchy-characteristic evolution, a number of systems are currently under inves- 
tigation. The motivation behind these new formulations of the initial value problem is to 
capture the advantages that each approach exhibits and at the same time avoid some of 
the corresponding drawbacks. The cIVP permits the construction of integration algorithms 
that allow null infinity to be included in a compactified grid, hence facilitating and clarify- 
ing the extraction of radiation. However, characteristic formulations generally break down 
in regions of complicated caustic structure, which are unavoidable in strongly asymmetric 
geometries such as those describing the merger of two black holes. Cauchy evolutions avoid 
this problem, yet they provide no natural way to impose conditions at the outer boundary. 
Since each method operates successfully precisely in the region where the other has its short- 
comings, an appropriate matching of the two initial value formulations promises an effective 
approach to the outer boundary condition and the caustic problem. 

The general idea behind the Cauchy-characteristic matching (CCM) is not entirely new. 
An early mathematical investigation exhibiting unions of space-like and characteristic sur- 



faces was given in ||12|| . Regarding general relativistic systems, a discussion of the potential 
of the method appears in The concept of a null exterior attached to a Cauchy evolution 
appears also in connection with perturbative approaches to the outer boundary problem, 



e.g., in [T9|] as well as in M. Yet, only recently P-[TI|J^], has the concept been carefully 
explored with respect to its practicability. A detailed study of the stability and accuracy 
of CCM for linear and non- linear wave equations has been presented in [f||§, illustrating 
its potential for a wide range of wave systems. The numerical investigation of cylindrically 
symmetric solutions of the Einstein equations has also been carried out [|IO| ,|il"|l . 

The objective of this paper is to develop and carefully calibrate the CCM method for 
asymptotically flat, spherically symmetric space-times, which are evolving in the presence 
of a minimally coupled, self-gravitating, massless scalar field. This is an initial step towards 
developing a general method applicable to the full Einstein equations. Research on this topic 
is stimulated and guided by the requirements of the Binary Black Hole Grand Challenge 
Alliance, a major collaboration aimed at the investigation of the merger of two in-spiraling 
black holes. The CCM approach will provide, in this context, both boundary conditions 
and radiation waveform extraction. The model problem investigated herein captures some 
essential aspects of the general system, including wave propagation on dynamical back- 
grounds and black-hole formation. The dynamics of the model is governed by the coupled 
Einstein-Klein-Gordon equations in spherical symmetry, 

G ab = K (V a $V b $ + Lg ab ) , (1) 

and 

V a V a $ = 0, (2) 
where $ is a scalar field whose Lagrangian is 

L = -^V a $ V a $ (3) 

and k = 8tt. (Units are chosen so that G = c = 1). Space-time indices are denoted with 
Latin letters (a, b, c, . . .) and space indices with (i, j, k, . . .). 
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This system has been studied intensively in the last five years both analytically || and 
numerically [0,0,|],|l3],[7| . The existent understanding of the system is quite appropriate, as 
the exploration and calibration of the combined evolution is done in a non-trivial geometric 
setting, yet with good basic knowledge of the expected physical behavior. 

The interior foliation and the associated integration algorithm are presented in Section [II]. 
The time integration is analogous to that used by Choptuik ||, the main difference being 
a different choice of boundary conditions. Section [JTJ describes the characteristic evolution 
algorithm. The method closely follows [16|, with minor changes to accommodate the match- 
ing interface. In Section [TV] the geometric concepts underlying the matching theory are 
discussed, and a number of simplifying assumptions are put forward, which lead to a simple 
set of matching conditions for the problem at hand. The implementation of the combined 
evolution, along with validation tests, and numerical experiments covering a wide range of 
methodologies and physical parameters are given in Section [V|. 



II. INTERIOR DOMAIN: CAUCHY EVOLUTION 

For the interior (Cauchy) domain M~ of the four- dimensional space-time, Einstein's 
equations will be written using a standard 3+1 ADM decomposition. Generally, a foliation 
of space-like slices (hypersurfaces) is constructed and labeled by a scalar function r. The 
unit normal to these slices is n a . By construction, n a = -«V a r, where a is the lapse 
function. The intrinsic metric on each slice is then given by 7 a ft = g a b + n a n . The vector 
an a connects the slices of the foliation; however, this time vector is not unique. In general, 
the vector t a can be chosen as t a = an a + f3 a where j3 a is the shift vector and j3 a n a = 0. 
Adopting this decomposition, the four dimensional metric element in the case of spherically 
symmetric space-times reduces to the 3+1 form 

ds 2 = (-a 2 + f3 2 e- 2A )dt 2 + 2/3 dt dp + e 2A dp 2 + p 2 e 2B (d6 2 + sin 2 6 d(J) 2 ), (4) 

where f3 = (3 p is the only non-zero component of the shift vector, and the metric on the 
space-like hypersurfaces takes the form 7^ = diag(e 2A , p 2 e 2B , p 2 e 2B sin 2 9). All the metric 
coefficients are functions only of p and t. The dynamical quantities in the 3 + 1 initial value 
formulation are the intrinsic metric 7^, and the extrinsic curvature K ao , of the slices. In 
spherical symmetry, the extrinsic curvature has only two independent components: K l - = 
di&g(H A ,H B ,H B ). 

The spatial projections of Einstein's equations lead to the evolution equations for the 
spatial tensors 7^' and Kj = j lk Kkj, 

A tt = -aH A + e- 2A ((3, p - (3A jP ), (5) 

B tt = -aH B + e~ 2A !3 Q + , (6) 

H A ,t = e~ 2A [aR pp + (3H Ap - a tPP + A if> a p ] + aH A K - n[e' 2A {^ >p )% (7) 

r / 2A-2B\ fi \" 

H Bjt = e~ 2A ai-—^\R ee + i3H B>p -a p i- + BA + aH B K, (8) 
where K denotes the trace K\ of the extrinsic curvature. 
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The time-time and time-space projections of the Einstein's equations yield respectively 
the Hamiltonian and momentum constraints 

-R pp + —Roe + e 2A H B (H B + 2H A ) = - K ($ p ) 2 + e 2A U 2 (9) 



and 



2 p 2 " x " 2 



H B , P + + (^5 - #a) = \^, P , (10) 
where i? pp and Rgg are the 3-Ricci components given by 

\r pp = -B >pp + -(A p - 2B, p ) + B, p (A p - B tP ), (11) 

p 2A-2B i i 

-Ree = ~B iPP + -(A p - AB >P ) + B jP (A >p - 2B jP ) + -(e 2A ~ 2B - 1). (12) 



Finally, the 3+1, first order in time, form of the Klein-Gordon equation (§) is 

$, f =«n + e- 2 > iP , (13) 



n,t = e" 2A |a ^ (p 2 $, P ) p + (2B ;P - A^ <p + e 2A KU 



+ a iP $ iP + PU P . (14) 



The geometry evolution equations (fj)-© together with the matter evolution equations 
(^3|) -(|l4]) constitute an initial value problem for the quantities A, B, Ha, H b , n. The 
initial data for the I VP must satisfy the constraints (p|)-(PI)|). Obtaining the development 
of the initial data requires, furthermore, the specification of the gauge functions a and (3. 
Ideally, the time integration of the IVP defined above should put no restrictions on either the 
form of the metric or the gauge functions. Yet the use of spherical coordinates, mandated by 
the Killing symmetries (and also computationally efficient), significantly limits the freedom 
to choose the space-time slicing. For example, the simplest possible gauge, the geodesic or 
synchronous gauge (a = 1,(3 = 0), leads to a coupled system of equations for the geometric 
variables A, B, Ha, H b . Formulating this problem as set of second order equations for A, B 
reveals a dynamical structure that is not that of a wave system. In practice, this leads to 
considerable difficulties in preserving the appropriate regularity of the geometry near the 
origin. Such regularity is analytically ensured by cancelations of 1/r terms both explicitly 
in the right-hand side of the evolution equations and implicitly by the enforcement of the 
constraints. 

A gauge choice that overcomes the regularity problems at the origin is the radial gauge, 
in which p is chosen such that the area of each dt = dp = sphere is equal to 47rp 2 , which 
leads to the condition B = 0. The integration procedure outlined below also assumes a 
vanishing shift condition, but can be generalized for arbitrary shift. 

The line element @ with B = /3 = becomes 

ds 2 = -a 2 dt 2 + e 2A dp 2 + p 2 {dB 2 + sin 2 6d<p 2 ). (15) 

As an immediate consequence of B = (3 = 0, equation @ implies that H B also vanishes, 
and equation (H) for H B t gives a condition on the lapse function, 
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a, P - 



A. p + -(e 2A -l) 
P 



a = . 



Furthermore, equations (|5]) and (|7|) reduce to 



A, = -aH A , 



-2A 



"Ap _ a ,PP + Ap a ,P 



+ aH A - Ke- 2A ($, P ) 



(16) 

(17) 
(18) 



With the conditions B = Hb = 0, the Hamiltonian and momentum constraints are now 
decoupled equations relating the metric function A and its associated extrinsic curvature 
H A to the scalar field energy density and current respectively: 



A. + 



,2A 



i) = p) 2 + e 2A n 2 



(19) 



Ha 



1 



(20) 



Finally, the scalar field equations, in this gauge, take the form 



$ t = all 

II t = ae~ J 



.2 y** 



A 

It* 



p $. + e A H A Il 



(21) 
(22) 



In summary, the lapse function is determined by the condition (|T^) . The remaining grav- 
itational variables (A and H A ) can be constructed from a subset of the evolution equations 
( |Hj - |T8| ) and constraint equations ( |T9"| - pO"|) . Different integration schemes can be designed de- 



pending on which two equations are chosen from (|T7| - |20|) to solve for A and H A . For instance, 
a free, unconstrained evolution consists of using equations (|I7]) and (|TE|). On the other hand, 
a fully constrained evolution would require solving the Hamiltonian constraint (|T9|), with the 
extrinsic curvature H A computed from the source terms using (pup. Alternatively, a mixed 
scheme can be followed, with H A again given by (|20|), while the metric variable is updated 
using (|17D , subject again to the lapse condition. Still other alternative schemes involving 
combinations of ( ]T7j - pO"|) are possible. A fully constrained evolution is adopted here since 
it has been used in earlier accurate calculations of scalar wave collapse |J and facilitates 
enforcing a regular metric boundary condition at the origin. 

In order to complete the IVP, specific boundary conditions for the scalar field variables 
$, II at both ends of the integration domain must be prescribed, along with integration 
constants for the hypersurface equations (|I9]) and (|16"D. The scalar field variables must be 
finite at p = 0, and thus the appropriate boundary conditions at this point are p$ = pU = 0. 
The outer boundary conditions for these variables are well understood only in the limit 
p —>■ oo, where rigorous outgoing wave conditions exist. Imposing a boundary condition at 
any finite p involves a certain physical approximation. Achieving a complete solution of the 
IVP without such additional assumptions is indeed the main focus of the CCM program, 
and the details are given in section [TV]. 
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The prescription of the integration constants for the hypersurface equations at some 
point po relates the labeling of the time and radial coordinates with the proper time and 
radial distance measurements of a privileged observer. Two obvious choices are geodesic 
observers, either at the center of symmetry, or at the outer boundary. While the choice does 
not have any consequences on physical observables, it is more natural for our integration 
procedure to assume a geodetic observer at the center of symmetry, and hence to require 
that A = a = there. In fact, this choice is well suited for the study of critical collapse as 
well, where the self-similar critical solution || manifests itself in terms of the proper time 
of an observer at the origin. 



For the exterior characteristic evolution, a family of null hypersurfaces u = const is 
introduced, emanating along the outward normals to the cross-sections of the matching 
world-tube and labeled by u = t. The outgoing null rays are parametrized by an area 
coordinate r, with r = R at the matching worldtube. The coordinate x = r/(R + r) is 
introduced for purposes of compactification, so that null infinity is located at x — 1. (In the 
numerical simulations to be presented later, we set the scale so that R = 1 and consequently 
x = 1/2 at the matching worldtube). In the null coordinate system, the line element in the 
exterior region has the form 



III. EXTERIOR DOMAIN: CHARACTERISTIC EVOLUTION 



ds 2 = -e 2X du(—du + 2dr) + r 2 {d0 2 + sin 2 9d(p 2 ), 



(23) 



r 



where the metric functions A and V depend only on u and r. 
The hypersurface equations for A and V are 



A, r = 2vrr($ r ) 2 



(24) 



and 



(25) 



The scalar wave equation in the characteristic region takes the form 



2(r$, tt ), r = -(V$, r ), r 



(26) 



In terms of the intrinsic metric of the (u, r) sub-manifold, 



r\ a pdx a dx^ = —e 2X du( — du + 2dr), 



(27) 



where the Greek indices take the values (0, 1), this reduces to 



where g = r$ and is the D'Alembertian associated with ^ Q/ g. 




(28) 
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The integration of the system (p4|)- (p6D proceeds with the specification of initial data 
$(wo,r) for r > R on the initial null cone Uq. The hypersurface equations are integrated 
radially and furnish compatible geometric functions A and V. This in turn allows the time 
integration of the scalar field equation, which provides new data in the neighborhood u^ + du 
and completes the integration cycle. The extra information needed to evolve the space-time 
is the characteristic initial data R) at the inner boundary, as well as the boundary 
values X(u,R) and V(u,R). A characteristic integration algorithm for (EH) may be based 



upon the null parallelogram made up of incoming and outgoing radial characteristics [17,15 



IV. CAUCHY-CHARACTERISTIC MATCHING 

The match is performed across a three-dimensional, time-like hypersurface E (world- 
tube), which divides the four-dimensional space-time into two disjoint sub-manifolds M + 
(Characteristic exterior) and M~ (Cauchy interior). Each sub- manifold is endowed with 
a metric g^ b that induces a unique intrinsic geometry at E. Independent coordinate charts 
{x } 1 * 1 are introduced in M^. Let s a be the space-like, unit normal to E directed from M~ 
to M + . The metric intrinsic to E is then given by 



lab 



g a b - s a s b , (29) 



where h b is the projection operator into the subspace E. The second fundamental form, or 
extrinsic curvature of E is defined by 

o ab = Kh d b w {cSd) . (30) 

If v a is the time-like, unit tangent to E, the metric h ab has the further decomposition 

Kb = -v a v b + R 2 q a b, (31) 

where q ab is the unit sphere metric. The metric tensors h ab , q ab and the vectors s a , v a satisfy 
the following orthogonality conditions: s a h ab = s a q ab = v a q ab = s a v a = 0. 

The metric continuity requirement can be recast as conditions on the induced norm of 
the tangent vector v a on E, 



v a v b g ab \- = [v a v b g ab ] + (32) 



and the surface-area radius, 



[q ab gab]- = [q ab gab} + . (33) 



In addition to the continuity of the metric across E, in order to prevent sheet discontinu- 
ities (singular hypersurfaces), the following conditions on the extrinsic curvature must be 
imposed: 

[v a v b e ab ]- = [v a v b e ab ] + (34) 

and 
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[q ab e ab ]- = [q ab eab} + - (35) 

Equations (|3^)- (|35|) should be understood as a limit process. Namely for any tensor A ab , 
the continuity condition [A a{) ]~ = [A ab ] + implies 

Q hm_ A~ b (Q) = o lim + A+ (O), (36) 

where G S 1 * 1 and Q — > P~,0 — > P + , through M~,M + , respectively. The continuity 
condition ( |36"D has meaning only if there exists a mapping Ai : E~ — > S + which transforms 
tensor components between the two coordinate systems. 

In the spherically symmetric case under consideration, the coordinate systems {t, p, 6, 4>} 
in M~ and {u, r, 6, <p} in M + are introduced. Because of the symmetry, the same coordinates 
9 and (f) are used everywhere throughout M^. This is not the case with the radial coordinate. 
In general, p and r are not the same. The line elements in are given by (|i~5"D and (p3|). 
respectively. The tangent vector to S in these coordinates is 

«» = -=i»=| ( * ,A0 ' 0) inM " (37) 
dr \ (it, r, 0, 0) in M+. 1 J 



In this section, a dot denotes derivative with respect to the proper time r along the world 

dr 

fae A (-p,i,0,0) inM 



tube E. That is / = $- — v a d a f. From w a s a = and s a s a = 1, it follows that 



2A. • • (38) 



e 2A (-r, it, 0,0) in M+, 



and 



8a= (a- 1 e A {p,a*e- A i,Q,0) inM" 
\ (-«, r + iiV/r, 0, 0) in M+, 

The continuity conditions on the metric, (^) and (|3"3"D, can now be rewritten as 

- a 2 i 2 + e 2A p 2 = -e 2X (ii 2 V/r + 2iir) (40) 

and 

p = r = R, (41) 

respectively, where R is the surface-area radius of the matching surface. R is taken to be 
constant on S, hence the matching surface is invariant under the orbits of the spherical 
symmetry. 

An additional assumption adopted here is that the matching surface does not move in 
coordinate space (f = p — 0) and that the coordinate time in M + and M~ agree in E 
{u = t). The metric continuity condition ( f40"D then reduces to 

a 2 = e 2 ^, (42) 
with all functions evaluated at the matching surface. 
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Similarly, the continuity conditions on the extrinsic curvature, ([34]) and (|35|), yield re- 
spectively 

[q ab V a s b }~ = [q ab V a s b } + (43) 

and 

[s a a a }~ = [s a a a } + , (44) 

where a a = M b V;,M a is the world tube "acceleration." In spherical symmetry, condition ( |43"|) 
reduces to 

[s a V a R} + = [s a V a R]-, (45) 
which can be rewritten using ( |3~9D explicitly in terms of the metric functions as 

T* = Ti> (46 » 

at the matching surface E. 

The condition (44) yields an equation for t from which the coordinate times, t or u, at 
E can be solved in terms of the proper time r. Alternatively, t(r) and u(t) can be obtained 
from the normalization condition of v a of s a . 

The matching of the scalar field variables across the two coordinate domains must ensure 
that neither the field values nor the field derivatives exhibit jump discontinuities on the 
interface. That is, 

[•]- = [*]+ (47) 

and 

[Z a V a $]- = [rV a $]+ (48) 

for any vector l a . It is convenient to choose l a = s a + v a , which is an outgoing null vector at 
E. In this case, the continuity equation at the matching surface becomes 

<9$ A d$ Vd$ 

m +ae ir P = R^- (49) 

In the limit R — * oo, for an asymptotically flat space-time with no incoming radiation, 
equation (|49|) reduces to the familiar Sommerfeld condition if the right-hand side of (|49D is 
set to zero. However, unlike the Sommerfeld condition which is only valid asymptotically, 
equation fl49|) is an exact relation, valid at any distance. 

For the momentum variable II, the procedure is different, to avoid imposing continuity 
on higher derivatives at E. Starting from the definition II = v a V a &, it is rewritten as 

n + s a V a $ = Z a V a $ (50) 

which leads to the condition 

dp dr 

The system of equations fl4~^ ) ,fl46| ) ,([49|) and (|5T|) completes the specification of the match- 
ing interface. 
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V. TESTS AND RESULTS 



A. Stability and accuracy tests 

The discretization algorithms in both domains, as well as the implementation of the 
continuity conditions are all constructed by replacing derivatives by second order accurate, 
centered finite differences. The matching surface lies at a fixed coordinate location, which 
is a grid point of both the p and r coordinate grids. This simple scheme leads to long 
term numerical stability, which is not an automatic feature of a matching algorithm [|J. 
Here, long term stability is defined as the bound evolution of initial data over time periods 
large compared to the light crossing time of the inner computational domain. This stability 
requirement is stronger than classical Von-Neumann stability, which requires bound local 
propagation of linearized modes. Such stability is becoming more and more important 
in numerical relativity as the desired integration times become longer. The CCM code 
developed in this work exhibited stability for any practical evolution time. 

In Fig. m a typical matched evolution is shown for initial data with a relatively small 
M/R ratio (0.08). The functional dependence of the initial data is given in this case by a 
Gaussian, 

r$ = \e~^-^l^ 2 (52) 

with A = 0.0225, r c = 2.0 and o = 0.1. The left column shows snapshots of the evolution 
in the Cauchy domain, with the radial coordinate p running from the origin p = to the 
matching radius p — 1. The right column shows the corresponding null evolution, with the 
compactified radial coordinate x running from the matching radius x = 0.5 to null infinity 
at x = 1. The field variables illustrated in all snapshots are p$ and r$, with the scale being 
constant for plots in the same row. 

The first row includes a few snapshots that follow an imploding pulse as it crosses the 
matching surface and propagates into the Cauchy domain. Note that at this instance, an 
inaccurate matching scheme would create back-reflection which would immediately register 
at null infinity. The apparent widening of the pulse as it enters the Cauchy region is due 
to the different radial functions used in the two coordinate systems as well as the doubling 
of the local propagation speed in the Cauchy sector. The second row demonstrates the 
propagation and reflection of the (marginally sub-critical) pulse off the origin. Strong non- 
linear distortions occur there, while the leading part of the pulse is already crossing the 
matching surface and radiating to null infinity. In the third row, the peak of the pulse 
is propagating outward across the matching surface; and the solution finally decays in the 
fourth row, as the trailing parts of the pulse cross S. The very slight curvature of the pulse 
in the null region, as the peak amplitude crosses the matching surface, indicates a small 
amount of backscattering occurring at this time. 

Second order convergence of all computed quantities to a limiting value is readily verified, 
e.g., by monitoring the final field configuration for a sequence of successively refined grids. 
A more physically intuitive test, the conservation of the total energy of the system, is 
a powerful probe into how well the discretization of the Einstein equations preserves the 
additional differential structure encoded in the Bianchi identities. A test of the absolute 
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convergence of the energy residual AM is performed for that purpose. The energy residual 
between two time levels Ui,U2 is defined as 

AM(ux,u 2 ) — M(ut) - M(u 2 ) + / P(u) du (53) 

J Ml 

where M(u) is the Bondi mass content of the space-time slice at time u while P{u) is the 
power flow at infinity AM is identically zero in the continuum limit and thus must converge 
to zero appropriately as a function of the discretization length-scale A. 

The explicit forms of the Bondi mass M and the radiated power P in terms of metric 
quantities are given by 

(54) 



M(u) = -e- 2 M-V 
2 r 



P{u) = -^e~ 2H Q 2 u (55) 

where H{u) = lim^oo A(m, r) and Q{u) = lim^oo r$(w, r). 

In Fig. [2] the convergence of the energy residual is demonstrated. A sequence of ap- 
proximate solutions with progressively finer resolution are obtained. The initial data (|52|) 
are prescribed in the null sector, while the Cauchy sector is taken to be flat initially. The 
Bondi mass is computed at the initial Bondi time u\ = and at a fixed final Bondi time 
U2 = 4.0, while the power integral is accumulated (to second order accuracy) at each integra- 
tion step. The second order convergence of the energy residual, consistent with the second 
order discretization of the component algorithms and the matching interface, is evidence of 
a successful matching of the two evolution schemes. The computational error in the mass 
and, most importantly, in the radiated power is directly controlled by the grid spacing A, 
and in fact diminishes as A 2 . 



B. Computing a space-time with two alternative foliations 

A comparison of the numerical solution obtained by two, considerably different foliations 



of our model spacetime is illustrated next. First, a calibrated characteristic code is used 
to obtain a global evolution of initial data corresponding to an incoming wave with support 
outside a radius X\ = r 1 /(l + r x ) (Fig. ^). This one-dimensional characteristic evolution 
will supply the waveform of the outgoing radiation coming back out to future null infinity. 
Since the initial data have compact support outside x±, the space-time portion delimited by 
the initial time (u = 0) surface, the origin world-line x = and the incoming null cone C 
(beginning at x\, will be flat. 

Next, the CCM code, with a matching radius at x = 1/2 < x\, is used to evolve the 
same characteristic data, along with flat-space Cauchy data on the initial time (t = 0) 
surface (Fig. |j). The incoming pulse enters the Cauchy region across the matching surface 
in the inward direction, then gets Cauchy evolved until it leaves the matching surface in the 
outward direction and ends up at future null infinity. This test compares the waveforms at 
null infinity produced by a global null code and by a CCM code. Theoretically, the general 
covariance of the equations guarantees that the output should be identical. In practice, 
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this test checks a combined algorithm (CCM) against a well calibrated scheme, i.e., the one 
produced by the global null algorithm. 

The waveforms are compared at null infinity, as would be measured by asymptotic 
(Bondi) observers. However, the time coordinates are considerably different for the two 
evolution schemes. In the CCM approach, the central time parametrizes the space-like foli- 
ation with the proper time of an observer at the center of symmetry. This parametrization 
subsequently also labels the outgoing null cone foliation of the outer region of the space-time, 
with the synchronization performed at the matching radius. In contrast, the time parameter 
of the null foliation follows the proper time of the central observer, directly labeling the out- 
going null cones emanating from the center of symmetry. Before comparing the two signals 
{gccM{t), 9n{u)) they must be reparametrized according to the asymptotic Bondi time. 

In a Bondi frame, geodesic observers would measure an asymptotically Lorentzian line 
element, which in the characteristic coordinate system is 

ds 2 = -du 2 - 2dudr + r 2 {dd 2 + sin 2 Odcf). (56) 

The Bondi time u is related to the coordinate time u of a general null cone foliation by the 
factor 

jjj __ 2H(u) 

du 

defined following equation (|55|) . 

In Fig. H the signals at null infinity and their difference are shown as a function of Bondi 
time. The initial (characteristic) data are 

$ = Kr 2 e~ {r - Tc)2/a \mkr (58) 

with the parameters for this plot being A = 6 x 10~ 4 , r c = 3.0, a = 0.6 and k = 10. 
This value of A is just below the threshold of black hole formation (which occurs at about 
A = 6.125 x 10~ 4 ) and leads to the strong distortion of the signal in the second half of the 
pulse. The signals obtained with the two codes overlayed as functions of Bondi time show 
little difference to graphical accuracy, a manifestation of physical covariance and algorithmic 
compliance. The grid sizes used for this run where 500 + 500 points for the CCM code, and 
1000 for the null code. The relative difference between the two signals for those resolutions 
is at the level of 0.1%. The maximum absolute value difference between the two signals over 
the total integration time provides a strong and physically interesting norm. This norm 
converges to zero with a measured rate of 1.99, consistent with the anticipated second order 
convergence. 

The investigation is now extended to strong field phenomena, with the study of initial 
data that end up in the formation of a black hole. In the Cauchy region, black hole formation 
is signaled when e 2A — > oo and the function 2m /r = 1 — e~ 2A — > 1 at some radius Rbh- 
The mass of the black hole is then Mbh — %Rbh- Alternatively, the behavior of the 
lapse function, or the scalar field values themselves can be used for similar estimates. In 
comparison, the formation of the black hole in the null region is signaled by the infinite 
red-shift between Bondi time and coordinate time for all r > Rbh, and the resulting decay 
of further radiative loss to infinity. The final value of the Bondi mass is then equal to Mbh- 
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(57) 



In practice, the evolution of the system cannot proceed accurately long after the black 
hole has formed. The use of a surface area coordinate r leads to a coordinate singularity 
inside the horizon H and the premise of the validity of the finite difference approximation 
breaks down. This is reflected, for example, in the breakdown of the convergence properties 
of the algorithm. Nevertheless, the matching procedure is sufficiently accurate to allow an 
investigation of black hole formation quite close to the matching surface. (The mass of the 
black hole, and consequently the Mbh / R ratio is controlled directly by the parametrization 
of the initial data.) 

Fig. shows the evolution of a strong Gaussian pulse, leading to a black hole that has 
a mass only slightly smaller than the mass scale set by the matching radius (R = 1). For 
demonstration purposes, the radial coordinate p is also compactified here, and the field 
functions are shown as if they were defined on a single grid. The highly non-linear pulse 
creates strongly back-scattered radiation as it propagates inwards, but not enough to avoid 
a black hole collapse. The discontinuity of the field derivatives at the location of the horizon 
is prominent in the late time snapshot. 

Fig. [7] displays the effects of black hole formation on the time development of the wave- 
form and of the Bondi mass. The initial data for those runs are given by ( |58"D with the 
amplitude A increased progressively from 5.50 x 10~ 4 to 6.25 x 10~ 4 in three steps. Up 
until Bondi time of about 7.0, an increase in the amplitude of the initial data leads to a 
near linear increase of the initial mass (lower plot) and to a corresponding red-shifting of 
the waveform. Afterwards, there is a bifurcation point, at which data that lead to a black 
hole space-time produce a signal severely distorted with respect to that of non-collapsing 
data. The Bondi mass either drops to zero for an asymptotically Minkowskian space-time, 
or coasts to a constant value that represents the mass of the newly formed black hole. 

VI. CONCLUSION 

The model problem of the self-gravitating scalar-field provides a convenient framework for 
the study of the CCM research program in a controlled situation. It has been demonstrated 
here that the matched evolution is essentially transparent to the presence of the interface. 
The targets of (1) clearly identified radiative quantities and (2) physically accurate boundary 
conditions are achieved in this case with minimal computational and developmental effort. 
Direct application of the continuity conditions at the matching interface leads to a stable 
and accurate mixed evolution algorithm for all relevant Mj R ratios. Although idealized, this 
model problem captures many essential physical aspects of more generic asymptotically flat 
space-times. It is expected that the geometrical approach and the algorithmic methodology 
initiated here are applicable and useful in more realistic problems, and work is currently 
under way exploring this possibility. 
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FIGURES 



FIG. 1. CCM evolution. The left column depicts evolution in the Cauchy domain, the right 
column shows the corresponding null evolution. The field variables illustrated are p& and r$ 
respectively. For each row, the left hand scale gives the amplitude in the two domains. 

FIG. 2. Convergence test for the energy residual AM. The grid size Ap refers to the Cauchy 
sector. The sequence of successive higher resolutions maintains a fixed ratio of the null sector grid 
size (Ar) to that of the Cauchy sector grid size (Ap). The convergence rate for the demonstrated 
sequence of grid resolutions is 1.89. 

FIG. 3. The foliation of an asymptotically flat spherically symmetric space-time by outgoing 
null cones emanating from the origin. The initial data have compact support outside x\ so that 
the inner region of the space-time (bounded above by C) is flat. The vertical and horizontal dotted 
lines give the location of the matching radius and the initial data surface of the comparison CCM 
run. 



FIG. 4. The foliation of an asymptotically flat spherically symmetric space-time with a com- 
bination of space- like and null hypersurfaces. Initial data for the comparison with the globally 
null code are given outside the matching surface. Flat Cauchy data in the interior complete the 
specification of a physical system. 



FIG. 5. Comparison of null and CCM evolutions in the case of sub-critical initial data. The 
upper plot demonstrates the nearly identical signal at null infinity. The solid line shows the signal 
output of the CCM code, while the circles indicate the corresponding null code output at the same 
Bondi time. The lower plot highlights the difference in the signals. The signal difference for a 
CCM grid of 500+500 points and a Null grid of 1000 points, is illustrated by the circle line. The 
diamond line shows the second order reduction of the difference when the resolution of all grids is 
doubled. 

FIG. 6. CCM evolution and black hole formation. The evolution of sufficiently strong Gaussian 
initial data leads to the formation of an apparent horizon, in this case at about p = 0.8 (x = 0.45). 
This is just marginally inside the matching radius at p = 1 {x = 0.5). The incoming Gaussian 
pulse (right side of the plot) propagates inwards (to the left), crosses the matching surface and 
collapses, forming a cusp-like profile. 

FIG. 7. The distortion of the signal at infinity for a space-time that develops a black hole. The 
upper plot shows overlapping signals with a progression of initial amplitudes. The dotted signal 
corresponds to black hole formation. The lower plot shows the evolution of the total Bondi mass. 
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